%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stationary optical Bloch equations for a Lambda system,
%  n = [22 33 12 21 13 31 23 32],
%
%	   --- |2> 
%	   / \			
%	  /  --- |3>
%	 /	  		  
%  --- |1> 
%
% Susceptibility spectrum obtained when ramping the pump laser,
%  fixed dressing laser.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Bloch
    figure(1);			clf;
    G12=              1;                                                    % decay of populations
    G23=              1;
    g13=  [0:.3:1]*1;
    o12=          1*G12;													% Rabi frequencies
    o23=	      1*G23;
    D12=          -G12*4;                                                   % detunings
    D23=[-1:.02:1]*G12*4;
    d23=	          0;

    for kk=1:length(g13)
      G13=	    g13(kk);
      for ll=1:length(D23)
            d12=     D23(ll);
            n  = Matrix_L(G12,G23,G13,o12,o23,d12,d23);
            nt(kk,ll)=real(n(1));
            n5r(kk,ll)=real(n(5));
            n5i(kk,ll)=imag(n(5));
        end
    end

    figure(1);			clf;
    subplot(331),plot(D23/G12,nt','LineWidth',1.2);
    set(gca,'xLim',[-4 4],'FontSize',16,'LineWidth',1.3)
    xlabel('\fontsize{16}\Delta_{23}/\Gamma_{12}');
    ylabel('\fontsize{16}{\it\rho}_{22}');
    text(-7,.15,'\fontsize{16}(a)')
    subplot(332),plot(D23/G12,n5r(1,:)',D23/G12,n5i(1,:)','LineWidth',1.2);
    set(gca,'xLim',[-4 4],'FontSize',16,'LineWidth',1.3)
    xlabel('\fontsize{16}\Delta_{23}/\Gamma_{12}');
    ylabel('\fontsize{16}{\it\rho}_{12}');
    text(-7,.5,'\fontsize{16}(b)')

% susceptibility
% 
%     E0 =8.854e-12;
%     hP =1e-34;
%     n  =1e18;                                                             % atomic density
%     k  =2*pi/689e-9;														% wavelength
%     dge=sqrt(3*pi*E0*hP*G23/k^3);											% dipole moment
%     chi=n*dge^2/E0/hP/o12;                                                % susceptibility factor
% 
%     subplot(224),plot(D12/G12,nt5*chi,'r','LineWidth',1.2);
%     set(gca,'FontSize',16,'LineWidth',1.3)
%     xlabel('\fontsize{16}\Delta_{ge}/\Gamma_{ge}');
%     ylabel('\fontsize{16}Im {\it\chi}');

    DarkResonance;
end

function DarkResonance
    G12=              1;                                                    % decay of populations
    G23=              1;
    G13=	        .001;
    o12=          1*G12;													% Rabi frequencies
    o23=	      1*G23;
    D12=[-1:.03:1]*G12*4;                                                   % detunings
    D23=[-1:.03:1]*G12*4;
    d23=	          0;

    for kk=1:length(D23)
      d23=	 D23(kk);
      for ll=1:length(D12)
            d12=     D12(ll);
            n  = Matrix_L(G12,G23,G13,o12,o23,d12,d23);
            nt(kk,ll)=real(n(1));
        end
    end

    subplot(333),mesh(D23,D12,nt');
    set(gca,'FontSize',16,'LineWidth',1.3)
    xlabel('\fontsize{16}\Delta_{23}/\Gamma_{23}');
    ylabel('\fontsize{16}\Delta_{12}/\Gamma_{12}');
    zlabel('\fontsize{16}{\it\rho_{22}}');
    text(-6,7,.18,'\fontsize{16}(c)')
    shading flat
    axis tight;
end

function n = Matrix_L(G12,G23,G13,o12,o23,d12,d23)
% trace-reduced Liouville matrix for a three-level cascade system

    I=          1i/2;
    g12=     G12+G23;                                                       % decay of coherences
    g13=         G13;
    g23= G12+G23+G13;
    d13=     d12-d23;
    L12=g12/2+1i*d12;
    L13=.05+g13/2+1i*d13;
    L23=g23/2+1i*d23;

    b=[0 0 -I*o12 I*o12 0 0 0 0]';
    
    M=[
       -G12-G23	0		-I*o12	I*o12	0		0		I*o23	-I*o23
        G23		-G13	0		0		0		0		-I*o23	I*o23
        -1i*o12	-I*o12	-L12	0		I*o23	0		0		0
        1i*o12	I*o12	0		-L12'	0		-I*o23	0		0
        0		0		I*o23	0		-L13	0		-I*o12	0
        0		0		0		-I*o23	0		-L13'	0		I*o12
        I*o23	-I*o23	0		0		-I*o12	0		-L23'	0
        -I*o23	I*o23	0		0		0		I*o12	0		-L23
    ];

    n = -M\b;
end
